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Abstract 

We examine the reversible adsorption of hard spheres on a random site surface 
in which the adsorption sites are uniformly and randomly distributed on a plane. 
Each site can be occupied by one solute provided that the nearest occupied site is 
at least one diameter away. We use a numerical method to obtain the adsorption 
isotherm, i.e. the number of adsorbed particles as a function of the bulk activity. The 
maximum coverage is obtained in the limit of infinite activity and is known exactly 
in the limits of low and high site density. An approximate theory for the adsorption 
isotherms, valid at low site density, is developed by using a cluster expansion of the 
grand canonical partition function. This requires as input the number of clusters of 
adsorption site of a given size. The theory is accurate for the entire range of activity 
as long as the site density is less than about 0.3 sites per particle area. We also 
discuss a connection between this model and the vertex cover problem. 

Key words: Adsorption; random site surface; cluster expansion; vertex cover 
PACS: 05.70.Np, 05.20.-g, 



1 Introduction 



The adsorption of proteins and colloids at a liquid-solid interface is a key step 
in many natural and industrial processes such as nitration, chromatography, 
protein purification, immunological assays, biosensors, biomineralization and 
biofouling. In many of these situations, the surface of the adsorbent is hetero- 
geneous. For example, in immunological assays one typically employs colloidal 
particles that have been coated with proteins (antibodies) to bind with anti- 
gens that may be present in the sample. Similarly, in affinity chromatography, 
the adsorbent is synthesized by immobilizing certain affinity ligands on porous 
silica, agarose or synthetic polymers. Many of these applications can benefit 
from a quantitative knowledge of the amount of solute that is adsorbed as a 
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function of the bulk concentration, ligand density and distribution and solute 
properties such as the size. 

The modeling approach required depends on the nature of the adsorption. 
When there is a finite desorption probability, characterized by a non-zero des- 
orption rate constant, an equilibrium between the bulk and adsorbed phases 
will be established; rapidly if the desorption rate is large and more slowly for 
small desorption rates. The properties of the equilibrium state, in particular 
the adsorption isotherm, depend on the bulk phase activity. If the desorption 
rate constant is very small on the experimental time scale the adsorption is 
effectively irreversible and a different approach is required. 



The statistical mechanics of reversible adsorption on heterogeneous surfaces, 
particularly for gases on solids, has a long history. Pioneering work was per- 
formed by Hill (y]) and Steele (0) and it is still an active area of research 
0). There is also a well-developed literature on irreversible adsorption on 
homogeneous surfaces (0;S0). Some models specifically address irreversible 
adsorption on non- uniform surfaces (@; 0; 11) where macromolecules are 
represented as hard spheres that bind irreversibly to adsorption sites. In the 
simplest of these, the Random Site Surface (RSS) fef), the sites are represented 
by randomly distributed points. Adamczyk et al extended the basic model 
to the situation where the adsorption sites have finite dimensions. The adsorp- 
tion of colloidal (H 3) and nanoparticles (14) has been interpreted with these 
models and a similar hard sphere model was used to rationalize the adsorption 
of proteins to hydrophobic sites on mixed self-assembled monolayers ()15l). 



Although macromolecules such as colloids and proteins have a tendency to 
adsorb irreversibly, this is not always the case and it is certainly useful to 
understand the equilibrium behavior. It this article we therefore present nu- 
merical and approximate analytical results for reversible adsorption of hard 
spheres on the RSS. For irreversible adsorption on this surface the task of de- 
veloping a theoretical description was greatly simplified by the existence of an 
exact mapping to an analogous process on a continuous (homogeneous) sur- 
face (0). For reversible adsorption on the RSS surface, however, there appears 
to be no similar mapping that would allow us to exploit the known behavior 
of hard spheres on continuous surfaces 



HQ. 



In addition to its application to adsorption, the RSS model is also interesting 
from another perspective. Weight and Hartmann obtained analytical results 
for the minimal vertex cover on a random graph by a mapping to a hard sphere 
lattice gas (jlH llih . A vertex cover of an undirected graph is a subset of the 
vertices of the graph which contains at least one of the two endpoints of each 
edge. In the vertex cover problem one seeks the minimal vertex cover or the 
vertex cover of minimum size of the graph. This is an NP-complete problem 
meaning that it is unlikely that there is an efficient algorithm to solve it. 
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Fig. 1. Sample configurations of adsorbed hard spheres on the random site surface. 
A sphere may bind, centered, to a site as long as the nearest occupied site is at 
least a away. Left configuration: a = 1, N s = 200,$ = 0.380; right configuration: 
a = 10, N s = 2000, 6 = 0.650. In both cases A = 10000. 

The connection to the adsorption model is made by associating a vertex with 
each adsorption site. An edge is present between any two vertices (or sites) 
if they are closer than the adsorbing particle diameter. The minimal vertex 
cover corresponds to densest particle packings. Weight and Hartmann obtained 
an analytical solution for the densest packing of hard spheres on random 
graph using a replica-symmetric approach. Although the random graph is 
related to the RSS model, it is not the same as the former does not have 
a physical structure. Specifically, in a random graph each possible edge is 
present with a given probability, c. In the random site surface, on the other 
hand, two adsorption sites that are neighbors of a given site are more likely to 
be neighbors of each other than two randomly selected sites. This effect can be 
quantified by the clustering coefficient, which is the average probability that 
two neighbors of a given vertex are also neighbors of one another. 

In section 2 we define the model and the simulation procedure. Section 3 
compares numerical results for the maximum coverage with various theoretical 
estimates. The structure of the random site surface is discussed in section 4, 
and section 5 presents an approximate theory for the adsorption isotherms 
that applies at low adsorption site density. 



2 Model and Simulation 



The adsorbent surface consists of N s point sites that have been uniformly and 
randomly distributed on a surface of area L 2 . These sites are frozen in place. 
The adsorbate molecules are represented by hard spheres of diameter a. In 
order to adsorb, a hard sphere must bind, centered, to one site. An adsorbed 
sphere may cover any number of sites but only occupies (or interacts with) 
the one at its center. Any site that lies within a distance a of an occupied site 
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is unavailable. See Figure 1. It is convenient to introduce a dimensionless site 
density 

corresponding to the average number of sites in an area equal to the projected 
area of a sphere. The coverage is defined as 
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where n is the number of adsorbed spheres. If the adsorption is irreversible it 
was shown in Ref (0) that the coverage at time t is equal to the coverage on 
a continuous surface at time r, 



9(t;a) = 9 c (r), (3) 

where r and t are related by 

r = a{l-e~ t/a ), (4) 

which has the required property that as a — > oo, r — > i. Although the time- 
dependent coverage on the continuous surface c (t), is itself not known exactly, 
semi-empirical equations that provide an accurate representation of numerical 
simulations of the model have been proposed (0). It is not clear, however, how 
to develop a similar mapping to apply to the reversible (equilibrium) case. 

In the reversible binding case there is an equilibrium between a bulk phase 
containing adsorbate at activity A, and the thermodynamic properties of the 
adsorbed phase can be obtained from the grand canonical partition function. 
For a particular realization of the adsorption site distribution this is: 



~(A,Ag= inc* -/<,•**",•) n** 1 ]. ( 5 ) 

rt;=0,l i>j i 

where A = exp(/5/i) is the activity with (3 = 1/ksT and \x is the chemical 
potential. Site i is occupied (vacant) if rij = 1(0) and = 1 if sites i and j 
are closer than a and zero otherwise. 

The number of adsorbed molecules can be computed directly from the parti- 
tion function: 

, fd\nZ\ 



as can the fluctuation in the adsorbed number: 

< n 



< (Snf >= A ( 8 -^\ (7) 
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It is generally not possible to evaluate this partition function analytically, 
even in one-dimension. However, if all the sites are isolated fij = 0,Vi,j the 
partition function becomes 

Z(\,N S ) = (1 + \) N °, (8) 
which yields the Langmuir form for the fraction of occupied sites, 

9 = a-^- (9) 

1 + A v ' 

and fluctuation 

< {8nf >= Ns-^^- (10) 



We have used a Gillespie type algorithm (|20t l21[ ) to obtain the adsorption 
isotherms of the RSS model numerically. First, Ns sites are distributed uni- 
formly and randomly in a periodic unit cell. A list of the neighbors, i.e. those 
that lie within a distance a, of each site is then constructed. At each step of 
the simulation, the total event rate is R = + n/A where is the fraction 
of sites that are available, and n is the number of occupied sites. A waiting 
time is generated from an exponential distribution using t = — ln(^)/i? where 
£ is a uniform random number. A second random number is used to decide 
the type of event. If £ < 0/(0 + n/A) a new sphere is added to the surface 
at an available site selected at random. Otherwise, one of the spheres on the 
surface is removed. For a given value of A, the adsorption/desorption process 
is allowed to continue until the system reaches a steady state and then the 
number of adsorbed spheres is averaged over a large number of events (typi- 
cally one million). In the steady state this procedure is equivalent to a grand 
canonical Monte Carlo simulation. 

All of the simulation results reported here were obtained using 4000 adsorption 
sites. Several tests were performed with larger and smaller systems, but there 
was no noticeable system size effect. Up to 20 different realizations of the 
adsorption sites were used. The isotherms were generated by performing a 
series of simulations with increasing activity in the interval A m i n < A < A max . 
The final configuration of a simulation at given activity served as the initial 
configuration for the following simulation at higher activity. For each value 
of the activity n tria i events were simulated. In order to allow the system to 
equilibrate at each new activity, only the second half of the ntrial events were 
used in the calculation of the ensemble averages. 

The activity schedule was determined from the formula: 

A » 

Ai = A min (-^)" A_1 , i = 0,...,n A -l (11) 
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Fig. 2. Adsorption Isotherms from simulation. Results are averages over 20 realiza- 
tions of 4000 sites, a = 1,2,4 top to bottom. 

For the adsorption isotherms at low activities we used A min = 0.1, A max = 100 
and ri\ = 20 and n trial = 50000. 

As long as the activity is not too large, or whenever a < 0.8, the simulations 
are non-problematic in that the steady state is reached rapidly and repro- 
ducibly. For larger values of alpha at high actvities, some care is needed to 
ensure that the system reaches equilibrium. The parameters n\ and n tria i were 
incresed until further changes had no influence of the isotherms. In most cases, 
n\ = 50 and ntriai — 5xl0 5 were sufficient. 

Alt houg h we did not employ the method in this study, parallel tempering 
( 22l ; 23f) would be an alternative, and possibly more efficient, method at large 
values of the activity. 



3 Maximum Coverage 



The maximum coverage is obtained in the limit A — > 00. As as the site density 
approaches zero, all sites become independent and the maximum coverage is 
a. In the limit a — > 00 the surface becomes smooth and an adsorbing sphere 
can occupy any position. The densest configuration in this case corresponds 
to a hexagonal close packing of spheres with a coverage of 6> hcp = 7r/(2\/3) = 



0.90687... fl2J). 



To obtain numerical values between these two limits for a given value of a 
we performed simulations using A m i n = 10, A max = 100000 and ri\ = 50. Some 
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sample isotherms, showing 1/9 versus 1/A, are shown in Figure 2. When 1 / A is 
sufficiently small, the curve approaches a plateau. These maximum values for 
different a are presented in Fig 3 together with various theories including the 
Langmuir approximation, the cluster theory to second order (see section 5) and 
the maximum coverage in an irreversible adsorption process (the "jamming 
limit") @). 

Weight and Hartmann (flit 0) used a statistical mechanics approach to the 
minimum vertex cover problem on a finite connectivity random graph. By 
using a replica-symmetric approach, they showed that the average maximum 
fraction of occupied graph vertices is 

_ 2W(c) + W(cf 

"max 2^ ' \ ) 

where W{x) is the Lambert- W function defined by W(x)e w ^ = x and c 
is the average connectivity of the random graph. From Eq 1 and Eq 2 it is 
easy to show that v = N/N s = 9 /a. Finally, in order to connect c to a we 
note that in a random graph, the distribution of connectivities is given by a 
Poisson distribution with mean c. Therefore, the probability that a vertex has 
a connectivity of zero is exp(— c). Similarly, for a random distribution of points 
in the plane at a density a the probability that a point has no neighbors (closer 
than a) is exp(— 4a) (c.f. Eq 15). Therefore, we take c = 4a. Substituting these 
results in Eq 12, the maximum coverage on a random graph is 

#max = V(4a) + V(4a) 2 . (13) 
4 o 



It is evident that the Langmuir approximation of independent sites is very 
poor, even for quite small values of a. As expected, due to the inherent ineffi- 
ciency of irreversible adsorption, the RSA coverage underestimates the maxi- 
mum coverage for all values of a. The remaining theories, all of which are equi- 
librium, consistently overestimate the maximum coverage. The cluster theory 
to second order provides a good estimate up to a ~ 0.3 



4 Structure of the random site surface 



The connectedness of the adsorption sites is a strong function of the dimen- 
sionless site density, as can be seen in Fig. 4. For a given distribution of sites, 
increasing this parameter corresponds to increasing the size of the adsorbing 
spheres. At small values of a, the surface consists of isolated clusters of sites. A 
cluster consists of n sites all of which are closer than a to at least one other site 
of the cluster. For small values of a the majority of sites are isolated. Above 
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Fig. 3. Maximum Coverage as a function of the inverse site density. Squares: numer- 
ical simulation; upper solid line: irreversible adsorption; dashed line: cluster theory 
to second order; dash-dotted line: random graph; dotted line: independent site as- 
sumption. 

a certain value of a c ~ 0.8 percolation occurs and almost all sites belong to a 
giant cluster (Jiih . 

In order to develop an approximate theory for the adsorption isotherms, we 
need to know the expected number of clusters of given type as a function of the 
site density. Analytical expressions for isolated sites and pairs can be obtained 
by making use of the following well-known result. Consider a region of area 
A that contains a finite number N of sites at given positions. The probability 
that none of the remaining N s — N points lies in this region (in the large N s 
limit) is 

Po = e~ p ° A . (14) 

where p s = N s /L 2 is the site density. A given site is isolated if there are no 
other sites within a circle of radius a centered on the given site. Therefore the 
expected number of isolated sites (per adsorption site) is 

x x = e~ 4a . (15) 



The number of pairs of sites that are separated by a distance between r and 
r + dr is ^2irrp s dr. Now imagine two circles of radius a centered on these 
sites. Let A 2 (r) represent the area of union of these two circles If any other 
site lies within this area, the two sites do not form a pair. The probability 
that no other site lies within this area is exp(— A 2 (r)p s ) where of diameter a 
separated by r. So the probability that two sites separated by r form a pair is 
p s iir exp(— A 2 {r)p s )dr . Integrating over r and introducing the dimensionless 
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Fig. 4. Structure of the random site surface for different values of a. Two sites are 
connected if thev are closer than a. 
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Fig. 5. Fraction of sites belonging to clusters of various sizes. From top to bottom, 
left hand side: isolated sites, pairs, triplets type a; triplets type b. The dashed line 
shows the fraction of sites not belonging to any of these types 
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variables a and y = r/a leads to the final expression for the number of pairs: 



x 2 =4a / exp A 2 (y) ydy, (16) 

where A 2 {y) = 2{ix — cos~ 1 (y/2)) —y^Jl — y 2 /4. Although this integral cannot 
be expressed in terms of elementary functions, it is straightforward to calcu- 
late numerically. A similar procedure could obviously be used to evaluate the 
expected number of higher order clusters, although the task rapidly becomes 
cumbersome. 

Numerical results for isolated sites, and clusters composed of two and three 
sites are presented in Fig 5. There were obtained by averaging over 100 inde- 
pendent distributions of 4000 sites each in a square cell with periodic boundary 
conditions. The theoretical expression, Eq 16 is indistinguishable from the nu- 
merical results. 

For a < 0.1 the surface consists almost entirely of isolated sites and pairs, 
while for a = 0.5 about 50% of the sites belong to clusters of cardinality four 
or greater. For convenience, we also fitted the numerical results using simple 
functions. For pairs the numerical data is very well fitted by the equation: 

x 2 = 2aexp(-5.578a), < a < 1. (17) 

There are two kinds of triplets. In type "3a" two sites are neighbors of a third 
site but not of each other, while in type "3b" each site is a neighbor of the 
other two. The number of each kind is well described by 

x 3a = 3.249a 2 exp(-7.537a), < a < 1 (18) 

and 

x 3b = 1.519a 2 exp(-6.250a), < a < 1. (19) 

As can be see in Fig 6, these fitting functions provide an excellent approxima- 
tion of the numerical results. 



5 Approximate theory for the adsorption isotherms at low site 
density 



Assuming that the adsorption surface consists of isolated clusters of sites, the 
partition function may be expressed as 



) N i(Z 2 ) N *(Z 3a ) N ^(~ 3b ) N ^..., (20) 



10 



Fig. 6. Fraction of pairs, and triplets of type 3a and 3b (top to bottom left hand 
side) as a function of the dimensionless site density from numerical simulation. The 
fitting functions Eqs 17,18, and 19 are shown as dashed lines. 
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Fig. 7. Clusters with grand canonical partition function. A solid line connects two 
sites that are closer than a. Triplets of type "3a" and "3b" are shown left to right 
in the second row. 

where iVj = XiN s is the number of clusters of type %. This factorization is 
possible because adsorption on a given cluster does not affect any of the others. 
The adsorption isotherm is, from Eq. 6, 



<n>=AiV 1 \ + XN 2 +.. 



(21) 



Figure 7 shows the first few possible clusters and the associated grand canon- 
ical ensemble partition function. A solid line connecting two sites indicates 
that they are closer than a particle diameter and therefore cannot be simulta- 
neously occupied. The partition function is evaluated by explicit summation 
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Fig. 8. Isotherms for a = 0.2,0.1, top to bottom. The solid lines show simulation 
results (an average of 10 distributions of 4000 sites). The dotted, dash-dotted and 
dashed lines show the cluster expansion to first (x%), second (X2) and third order 
(x 3a + x 3b ), respectively. 

over all sites. As an example, for the first of the two triplets, 



S 3 a = E E Ea-mr^Xl-m^A" 1 

ni=0 n,2=0 ri3=0 



+n 2 +n 3 



1 + 3A + A 2 . 



(22) 



By substituting these results into Eq 21 we obtain 



9 = a 



A 

TTa 



X\ + 



2A 



2A 



x 2 + 



3A + 2A 2 
1 + 3A + A 2 



X 3a + 



3A 
1 + 3A 



x 3b + ... 



(23) 



We show the predictions of this equation to different orders compared with the 
simulation results in Fig 8. Successive approximations approach the simulated 
isotherm from below. Including clusters up to triplets described the numerical 
simulation data accurately for a = 0.1, but is already quite poor for a = 0.2. 
This rapid breakdown can be understood by realizing that truncation of the 
expansion means that one is failing to include all sites. 

An alternative approach is to assume that all sites belong to clusters up to a 
given order. Thus, if all sites are assumed to be isolated we have 



which gives for the coverage the Langmuir expression: 



N ' - (l + Xf% (24) 



e ' = tt TTX (25 » 
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Fig. 9. Isotherms for a = 0.3,0.2,0.1, top to bottom. Solid and dashed lines show 
the simulation results and the Langmuir theory, respectively. 

For A — ► oo this gives 9\ = a, i.e. all sites are occupied. When it is assumed 
that only isolated sites and sites with one neighbor are present, the number 
of the latter is estimated as N 2 = (N s - Ni)/2 = N s (l - exp(-4a))/2. The 
partition function is 



" — "1 —2 



(l + A) JVl (l + 2A) 



A J 2 



(26) 



which leads to a coverage of 

iViA 2iV 2 A 



9 9 = ^ 



N s \l + X 1 + 2A 



aA- 



f A(l + exp(-4a)) 
(1 + 2A)(1 + A) 



(27) 



The maximum coverage at this level of approximation, obtained by taking the 
limit A — > oo, is 



# 2 ,max = -(1 + exp(-4a)) = a - 2a 2 + 0(a 3 ). 



(2* 



At the triplet level, we estimate the total number of triplets as = (N s — 
Ni — 2N 2 )/3. Of these, a fraction y 3a = x 3a / {x 3a + ^3fe) are of type a and the 
remainder are of type b. This yields: 



a 



(y, 



36 



3A 
1 + 3A 



V3a 



3A + 2A 2 1 - x x - 2x 2 2A 
1 + 3A + A 2 ^ 3 + 1 + 2A 



X-2 



A 

1 + A 



(29) 
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Fig. 10. Isotherms for a = 0.5,0.3,0.2,0.1, top-to-bottom. The solid lines show 
the simulation results, while the dashed and dash-dot lines show the approximate 
theories 82 and #3, respectively 

The maximum coverage is then: 

#3,max = ac[(y 3a + 1)(1 - x x - 2x 2 )/3 + x 2 + xi]. (30) 

Figures 9 and 10 show the simulated isotherms compared with the alternative 
theory. Unlike the theoretical estimates shown in Fig 8, the estimates generally 
lie above the simulation results. It is clear that the Langmuir model (all sites 
assumed to be isolated) provides a very poor description, even at quite low 
site densities. There is not much difference between the second and third order 
cluster expansions. Both provide a satisfactory description of the isotherm for 
a < 0.3. 



6 Conclusion 

We have obtained numerical results for the isotherms of hard spheres adsorbing 
on a random site surface. For a given configuration of sites the maximum 
coverage is obtained by taking the limit of infinite bulk phase activity. We 
developed an approximate theory based on a cluster expansion of the partition 
function that is accurate up to moderate site density. 

Weight and Hartmann obtained an analytic solution for a related model of 
hard spheres on a random graph using replica theory. Their solution applies 
only to the densest particle packing. Generalization of their approach to finite 
activity or to the RSS surface does not appear to be straightforward. Perhaps 
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an alternative approach will allow a description of the thermodynamics for a 
wider range of site densities than is possible using the theory proposed here. 
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